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ABSTRACT 



We present the exact calculus of the gravitational potential and acceleration along the symmetry axis of a plane, homogeneous, polar 
cell as a function of mean radius a, radial extension Aa, and opening angle A<f>. Accurate approximations are derived in the limit 
of high numerical resolution at the geometrical mean (a) of the inner and outer radii (a key-position in current FFT-based Poisson 
solvers). Our results are the full extension of the approximate formula given in the textbook of Binney & Tremaine to all resolutions. 
We also clarify definitely the question about the existence (or not) of self-forces in polar cells. We find that there is always a self-force 
at radius (a) except if the shape factor p = aAip/Aa — > 3.531, asymptotically. Such cells are therefore well suited to build a polar 
mesh for high resolution simulations of self-gravitating media in two dimensions. A by-product of this study is a newly discovered 
indefinite integral involving complete elliptic integral of the first kind over modulus. 
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1. Introduction 

The structure and dynamical evolution of astrophysical discs 
is mainly governed by gravity. In this context, the numerical 
computation of the gravitational potential and forces of discs 
is of fundamental importance and the improvement in accu- 
racy, resolu tion and computing time remains an interesting chal- 
lenge (e.g..lMueller & Steinmetzl 1 1995: Matsumoto & Hanawa, 
2003: lLondrilldl2004tlHurdl2005tlJuselius & Sundholml 120071 



Li et all 12008k In many (if not all) hy drodynamical simu 
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lations of flat, self-gravitating discs (e.g.. [Huber & Pfenniger. 
120011; IZhang et all 120081: fearuteau & Massed l2008t) . the poten- 



tial is derived ever ywhere in the disc by mea ns of Fast Fourier 
Transforms (e.g.. iBinnev & Tremaind 119871) . This is an effi- 
cient techni que, which can be extende d to deduce directly ac- 
celerations dBaruteau & Masseil [2008). A weak point however 
is in the determination of the self-gravitating component (i.e., 
the effect of a cell on itself), which involves an improper in- 
tegral. An approximation o f this self-potential was proposed in 
IBinnev & Tremaina dl987l) . but for a cell in which the surface 
density va ries with the polar radius R as R~ 3 ^ 2 . In the same 
conditions, Baruteau & Masset (2008) argued that the radial ac- 
celeration is expected to vanish at the center of each computa- 
tional cell. This assertion is certainly correct asymptotically for 
homogeneous cells as their size becomes smaller and smaller. 
We stress that because Newton's law goes like R 2 , the self- 
gravitating component generally provides a significant contribu- 
tion to the total potential or force, and must therefore be treated 
as properly as possible, especially at high resolution. Any error, 
even relatively small, may introduce artifacts or biases in mod- 
els. 

In this letter, we report the exact expressions for both the po- 
tential and radial acceleration due to a homogeneous polar cell 
inside the cell itself and study the high-resolution limit. We re- 
call in Sect.|2]the integral expression for the potential along the 



axis of a polar cell which can be decomposed into two parts: a 
potential created by an homogeneous disc and the other a poten- 
tial due to horseshoe-like surface. Sections [3] and [4] are devoted 
to the determinatons of these two contibutions. The results for 
the polar cell (potential and acceleration) are obtained in Sect. [5] 
We discuss the case of high resolutions in Sect. [6] and show that 
there is always a self-acceleration except if the cell has a special 
shape. These results are important since they contribute to the 
reliability and accuracy of current simulations of self-gravitating 
media on a polar mesh. 

2. The polar mesh 

We consider a planar, homogeneous, polar cell with inner polar 
radius a\, outer radius ai = a\ + Aa, lower azimuth <p\, and upper 
azimuth (f>' 2 = <p\ + A<p as shown in Fig.[T] At any position (R, </>) 
in the plane of this cell, the gravitational potential is given by 
Newton's generalized formula: 



<KR, 



ri r4> 



adadip' 



Va 2 + R 2 - 2aR cos((f>' - . 



(1) 



where So is the surface density. This cell has a single axis of 
symmetry defined by <f> — \{<t>\ + 4>' 2 )- If we introduce the quan- 
tity: 

2V^ ... 

m = (2) 

a + R 

as well as the new variable ff such that 29' — n - ((f)' - <p), the 
potential along its axis of symmetry is given by: 



ifr(R) = -2GL 



dad6' 



(3) 



Vl — m 2 sin 2 6' 

where 49 — 2n - A(f>. The integral over 0' can be decomposed 
into two integrals with as lower bound. We then have: 
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da I — m\K.(m) 
\R 



F(9, m)], 



(4) 
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and 




geometrical mean <a> 



mean radius a 



Fig. 1. A polar cell with inner radius a\, outer radius «2, radial 
extension Aa = a-i - a\, and opening angle A<p. 



where 



F(8, m) 



-f 

Jo 



dff 



vT 



(5) 



m 2 sin' 0' 



is the incomplete elliptic integral of the first kind, m is the mod- 
ulus, 6 is the amplitude, and K(m) = F(^,m) is the complete 
elliptic integral of the first kind. The interpretation of this rela- 
tion is presented in Fig. [2] and follows from the superposition 
principle: the first part of the integral of Eq. [4] is the potential 
due to a flat disc, namely: 



C, F 

i da —n 



M R ) = -2GS I da^\-mK{m), (6) 
and the second part is that of a flat, horseshoe-like surface: 



M(R) = -2GZo 



J da A /— i 



mF(8, m). 



(7) 




Fig. 2. A polar cell is the superposition of a disc and a horseshoe 
with negative surface density. 

We can easily derived a as a function of m from Eq.|2l as well 
as the corresponding derivative da/dm. We have, respectively: 



a I 1 + m! 
R 



m 



(8) 



, da 
2R 



(1 ±m') 



'\2 



■dm, 



(9) 



where m' = Vl - m 2 is the complementary modulus, +m' is 
for R < a (i.e., matter is located right to the position R) and —m' 
is for,R > a (i.e., matter is located to the left). From Eqs.© and 
(O, we deduce that the self-gravitating potential in the cell along 
its symmetry axis is: 



4r(R) = 4GZ {) R 



r mj - (1 +m') 3 
Jm, +m'm 3 



[K(m) - F(6, m)] dm. (10) 



We note that this expression is rarely seen in the literature. 

3. The homogeneous disc 

3.1. Classical derivation 

To calculate fid, one classically changes^ the modulus m of K 

by setting m - where x < 1. The new modulus x is either 
a/R = u < 1 or R/a = v < 1 depending on R. So, the potential 
due to an homogeneous disc is given by: 



MR) 

-4GE fl 



-£ 2 K(v)f, in region (I), 
JT uK(u)du — j* 2 ^jfi-dv, 

in region (II), 
J"" uK(u)du, in region (III), 



(12) 



where region (I) is for R < a.\, region (II) is for a\ < R < a-i, 
and region (III) is for R > «2 (see Fig. |2). Since the indefinite 
integrals in Eqs. (fT2l are knowrQ, a close-form expression for ip& 
can finally be deduced. It is: 



MR) 

-4GY R 



m: tad). 

[E(m) - u' 2 K(u)] l ui + [^fl 1 in (II), 
[E(m) - m' 2 K(m)1" 2 in (III). 

VL J mi 



(14) 



3.2. A new indefinite integral ? 

From a numerical point of view, it is preferable to use u and v 
as modulus of the complete elliptic integrals rather than m (see 
Sect. [6J. However, an interesting issue is raised if we perform a 
change of modulus in Eqs. (fT4l) . to restore m as the modulus of E 
and K. We find: 



MR) 

-4GE fl 



,r E ( m )WK(,n) -p in (I), 

1 i-m' ], r 



' E(m)-m'K(m) 
1+m' 
E(m)-m'K(m) 
l+m' 



l 1 + r gW^w r i n (ii), 

-im-> 

in (III). 

Jmi 

The transformation is (e.g. lGradshtevn & Rvzhikl . ll965h : 

'2^ 



(15) 



K 



1+x 



(\+x)K(x). 



(11) 



2 In particular, we have for any k < 1 (Gradshtevn & Rvzhik, [i965l) : 

J^K(/t)^ = -^ and J K(k)kdk = E(k) - k' 2 K(k), (13) 

where E is the complete elliptic integral of the second kind and 
k' = Vl - k 2 is the complementary modulus. 
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If we now compare these expressions with the first part of Eq. 
[TUl we conclude that: 

(1±k ' )3 K f^ _ E(k)±k'K(k) 

-K(k)dk — h , (16) 



P 



k 3 k' 



1 +k' 



for any modulus k and compl ementary mod ulus k! . This indefi- 
nite integral is probably new (lJeffrevLl2009l) . 

4. The homogeneous horseshoe 

To calculate if/^, we have attempted to proceed as for the disc, 
by applying a change of modulus. We failed, mainly because, 
as mentioned, the new amplitude a of F resulting from Eq. QT| 
now becomes a function of the modulus. This introduces a severe 
mathematical difficulty that we were unable to circumvent. The 
answer is however obtained by considering Eq. [T6lwith incom- 
plete elliptic integrals. Using the partial derivativqfl of F and E 
with respect to their modulus (keeping the amplitude constant), 
we find that: 



d_ 

dk 
where 



E(6, k) + k'F(8, k) 



1 +k' 



(i±*0 
k 3 k' 



■F(8,k)+H(0,k) (18) 



H(8, k) 



(l±g) 
2kk'A 



sin 28 



(19) 



and A = Vl - k 2 sin 2 9. In fact, J H(8, k)dk can be expressed 
exactly in terms of basic functions. Actually, after some algebra, 
we have: 



H(8, k)dk ■■ 



sin 28 



In- 



A + k' 



k! 1 
± In ■ 



1 + A 



(20) 



It follows from Eqs.([T0l>, (HUl, (O, and d20j> that the potential 
along the axis of the homogeneous horseshoe is given by: 



-4GE 7? 



(\ E(8 2 ,m)+m' F(8,m) 
1— m' 

in (I), 

E(8,m)-m' F(8,m) 



sin 29 
4 



In 



l+m' 



+ 



E(8,m)+m' F(8,m) _ sin2fl i (A+m')(l+A) 
(A-m')(l-A) 



1 -in' 

' in (ID, 

E(8,m)-m'F(8,m) 
l+m' 



(A+m')(l+A) ] m 2 
(A-m')(l-A)J mi ' 

l'"2 
Jl 



sin 28 i (A+m'Xl-A) ] 1 
4 (A-m')(l+A)J m 



(21) 



sin 28 i n (A+m')(l-A) 1 m ? 
4 (A-m')(l+A)Jm, 



in (III). 

As for Eq.[14j the expression for region (II) can be simplified. 

5. Potential and acceleration in the polar cell 

The potential inside a polar mesh and along its symmetry axis is 
found from Eqs.(fT5b) and d2"Tb). It is: 

if,{R) _ E(m 2 ) - E(8, mi) + m' 2 [K(m 2 ) - F(8, m 2 )] 

-4GT. () R ~ 1 - m' 2 

E(mi) - E(8, m x ) - m\ [K(mi) - F(8, mi)] 



1 + m'j 

sin 26» (Ai - m[)(l + Ai)(A 2 - m' 2 )(l - A 2 ) 
~T~ "(Aj +m' 1 )(l-A 1 )(A 2 + m^)(l+A 2 )' 



(22) 



In particular, we have (e.g. Gradshteyn & Rvzhik. [l965h : 



dE 
dk 



E-F 



and 



dF E- k' z F 



fcsin^cos ( 



k dk kk' 2 k ,2 VT^fcWe 

where E(9, k) is the incomplete elliptic integral of the second kind. 



(17) 




Radius R 



Fig. 3. The potential along the symmetry axis of a homogeneous 
polar cell with geometrical parameters a\ = 0.95, Aa = 0.1 and 
A0 = 0.1 (a = 1 and (a) m 0.99875). The minimum lies at 
^min ~ 0.9949. Also shown are the potential of the associated 
disc and horseshoe. 

where A 2 = 1 - m 2 sin 2 8 with i = {1,2}. This expression is 
exact. Figure [3] displays if/^, i/r^ and their difference if/ versus R 
for a typical cell. By deriving this formula with respect to R, we 
obtain the relation: 



dR R 



(1 ± m'Y 



{K(m) - F(8, m)} 



(23) 



which is the opposite of the gravitational acceleration gg along 
the cell axis. If we set m = Rja^ and tfr = ifr/ifr c where 

<A C = -2nG^a 2 [l - |l - (24) 

is the potential due to the cell at the origin R — 0, and 
\il^l{K(m)-F(8,m)}] m2 

t m Jmi 



S e (m) 



then Eq.l23l takes the form: 



(l- a) (,r -20 



diff iff 

— = —+Sg(m). 
dm m 



(25) 



(26) 



This is the general ization of the Or d inary Differential Equation 
(ODE) derived by iHure & Hersanj d2007l) to a polar cell with 
opening angle A<p = 2n - 48 (which becomes a disc for 8 = 0). 
This ODE can be easily solved numerically using standard 
schemes since if/ is known both at the center (see Eq. l24l i and 
at infinity where if/ = 0. The full generalization to polar cells 
where the surface density S is a power-law of the radius, i.e., 
S(a) = Zo(a/«2) s , seems straightforward. 

6. Approximations at high numerical resolutions 
and the "magic" shape factor. 

At high resolution, polar cells tends to become segments, 
squares, or arcs depending on the shape factor p of the cell de- 
fined by: 

dAcp = pAa, (27) 
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where a = i(ai + 02) is the center of physical cells. When both 
A<p <k 1 (i.e. 9 — > |~) and Aa <k a, w e can expand F(6, m), 
E(8,m), K(m), and E(m) over m' (e.g., ICarlson & Gusta fson. 
1985; Gradshteyn & RyzhiE 1965) to obtain an approximation 
for i// and gf> as a function of 7? and the cell geometrical param- 
eters (a, Aa, A<p). A radius plays a key role in current FFT-based 
Poisson solvers (see references in Sect. [TJ: this is the geometri- 
cal mean of the inner and outer radii, namely (a) - y]a\ci2- This 
radius is the center of computational cells when a radial logarith- 
mic mapping is applied to the entire physical mesh. For R = (a), 
the two modulus nt\ and mi are equal (and so are the comple 
mentary), and we have m' { - m' 2 
algebra, we find: 



4^ to first order. After some 

4a 



i/f((a)) ~ -2GE Aa 



arcsinh p + p arcsinh - 
P 



and 



gR((a)) ~ G^o~r- |2p arcsinh arcsinh p 

a \ p 



(28) 



(29) 



Figure H] compares exact values of if/((a)) and git((a)) given 
by Eqs.d22l and (1231 with approximate values obtained from 
Eqs.d28l and (l29l respectively, for three factors p - {0.5, 1,2}. 
We note that, for Aa/a < 10~ 3 , relative errors increase. This is 
due to the "low dynamic" of the variable m as it approaches unity 
(and Aa — > 0). In this case, it is advisable to use u and v as the 
moduli of elliptic integrals instead of m to restore the accuracy 
in Eq.l22l 
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Fig. 4. Top panels: exact potential and acceleration in a polar cell 
at radius R = (a) for three shape factors p (circles) as well as 
their approximations (lines) by Eqs.d28li and d29l . Bottom pan- 
els: decimal logarithm of the relative error. Also displayed are 
the results obtained with the "magic" value p mag ; c « 3.531. 



The potential always has its minimum inside the cell at a cer- 
tain radius R m [ n reagardless of the opening angle A(f>. This radius 
is found by equating the right-hand side of Eq.[23]to zero. It de- 
pends on the parameters a, Aa, and p. We expect R m [ n < a for 
large opening angles (or low azimuthal resolutions) and R mm > a 
for small ones (high azimuthal resolutions). When Ac/) — > 0, 
^min — > a as the polar cell becomes a rectangle. Since a + (a), 
there is in general a self-force at R — {a) at high azimuthal 



resolution. However, the gravitational acceleration exception- 
ally cancels out if the shape factor has the remarkable value 
Pmagic ~ 3.530937746935705. This is the single root of Eq.|29] 
Forp = pmagic we have tp({a)) * -5.922555527732666 GZ Aa. 
Thus, if the computational mesh (Aa, A<p) contains cells all with 
shape factor p = p mag ic, then there is no need to compute the 
gravitational force at R — (a), since it is zero in the high resolu- 
tion limit. 

More generally, since the two moduli m\ and mi appearing 
in Eq.|23]are functions of the ratio 02/01, it is always possible to 
find a shape factor for which the acceleration vanishes at a fixed 
value of R/ai. It can therefore be advantageous to build the polar 
mesh accordingly. For instance, let Nr and be the numbers 
of cells in the radial and azimuthal directions, respectively, and 
Ri a = Ri be the radius of the disc inner edge. Once the shape 
factor p has been numerically determined such that gu - (this 
is Pmagic in the high resolution limit), the radius R out = Rn„+i of 
the outer edge is found from the relation: 



jV>p + 7r 
N<pp - n 



N, t 



K in x|l + 2E£) (30) 
P AW 



With such a mesh, there is no gravitational acceleration at the 
center of computational cells. 

7. Concluding remarks 

We have reported exact expressions for the self-gravitating po- 
tential and acceleration in homogeneous polar cells whatever 
their shape. The high resolution limit has been analyzed and re- 
liable apj£02chTiatiomJ^e bera derived. We confirm the valid- 
ity of [Binnev & Tremaine ( 1987)'s prescription for the singular 
term ^(0, 0) to first order. We have shown that there is always a 
self-acceleration at the center of computational cells. However, 
this acceleration can cancel out at any radius by an appropriate 
choice of cell shape. This worifl contributes to the treatment of 
self-gravity in hydrodynamical simulations, and especially for 
discs. 
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4 A Fortran 90 package called PolarCELL is available at the follow- 
ing address: 

www . obs . u-bordeauxl . f r/radio/ JMHure/ intro2polarcell . html 



